library(specr)
library(fixest)
library(readr)
library(broom.mixed)
library(texreg)
library(dplyr)
library(fixest)
library(patchwork)
library(ggplot2)

# Cross sectional linear FE models that use aggregate data (2000-2019) for linguistic cleavages


source("funs_and_constants.R")

epro_year_group_clusters <- read_csv("data/epro_data_h2.csv")

mean_h2 <- epro_year_group_clusters %>% 
  filter(year >= 2000) %>%# only obs after 2000 because of EPR 2.0 coding
  group_by(gwgroupid, countries_gwid) %>%
  summarize(
    across(c(share_disagreement, disagreement, n_orgs, n_ed_langs, hhi_lang, linguistic_segments_bin, regaut,
             groupsize, incidence_flag, status_pwrrank, any_multiethnic, nightlight_total_pc_log, n_tek_groups),
           ~ mean(.x, na.rm = T)))

full_model_h2 <- feols(disagreement ~ linguistic_segments_bin + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + nightlight_total_pc_log + n_tek_groups | countries_gwid,
                       cluster = ~countries_gwid, data = mean_h2)

# n rels as IV
full_model_h2_n_langs <- feols(disagreement ~ n_ed_langs + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + nightlight_total_pc_log + n_tek_groups | countries_gwid,
                              cluster = ~countries_gwid, data = mean_h2)


full_model_h2_herf <- feols(disagreement ~ hhi_lang + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + nightlight_total_pc_log + n_tek_groups | countries_gwid,
                            cluster = ~countries_gwid, data = mean_h2)


# full_model_h2_no1clusters <- gfeolser(disagreement ~ linguistic_segments_bin + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + nightlight_total_pc_log + n_tek_groups | countries_gwid
#                        +(1+ nightlight_total_pc_log + n_tek_groups |gwgroupid) + (1+ nightlight_total_pc_log + n_tek_groups |countries_gwid), cluster = ~countries_gwid, data = epro_year_group_clusters %>% filter(n_orgs > 1), family = "binomial")

full_model_h2_max <- feols(share_disagreement ~ linguistic_segments_bin + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + nightlight_total_pc_log + n_tek_groups | countries_gwid,
                           cluster = ~countries_gwid, data = mean_h2)

full_model_h2_max_n_langs <- feols(share_disagreement ~ n_ed_langs + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + nightlight_total_pc_log + n_tek_groups | countries_gwid,
                                  cluster = ~countries_gwid, data = mean_h2)

full_model_h2_max_herf <- feols(share_disagreement ~ hhi_lang + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + nightlight_total_pc_log + n_tek_groups | countries_gwid,
                                cluster = ~countries_gwid, data = mean_h2)

screenreg(list(full_model_h2, full_model_h2_n_langs, full_model_h2_herf, full_model_h2_max, full_model_h2_max_n_langs, full_model_h2_max_herf),
          custom.header = list("Disagreement (1/0)" = 1:3, "Prop. Disagreement" = 4:6))


####
## 2-way FE models (lpm)
####


full_model_h2 <- feols(disagreement ~ linguistic_segments_bin + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
                       + n_tek_groups | year +countries_gwid, cluster = ~countries_gwid, data = epro_year_group_clusters)

# n rels as IV
full_model_h2_n_langs <- feols(disagreement ~ n_ed_langs + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
                              + n_tek_groups | year + countries_gwid, cluster = ~countries_gwid, data = epro_year_group_clusters)


full_model_h2_herf <- feols(disagreement ~ hhi_lang + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
                            + n_tek_groups | year + countries_gwid, cluster = ~countries_gwid, data = epro_year_group_clusters)

full_model_h2_max <- feols(share_disagreement ~ linguistic_segments_bin + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic  + n_tek_groups | year + countries_gwid,
                           cluster = ~countries_gwid, data = epro_year_group_clusters)

full_model_h2_max_n_langs <- feols(share_disagreement ~ n_ed_langs + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic  + n_tek_groups | year + countries_gwid,
                                  cluster = ~countries_gwid, data = epro_year_group_clusters)

full_model_h2_max_herf <- feols(share_disagreement ~ hhi_lang + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic  + n_tek_groups | year + countries_gwid,
                                cluster = ~countries_gwid, data = epro_year_group_clusters)


screenreg(list(full_model_h2, full_model_h2_n_langs, full_model_h2_herf, full_model_h2_max, full_model_h2_max_n_langs, full_model_h2_max_herf),
          stars = stars)

# removing colinear models
texreg(list(full_model_h2_n_langs, full_model_h2_herf, full_model_h2_max_n_langs, full_model_h2_max_herf),
       stars = stars,
       custom.coef.map = coef_map,
       custom.header = list("Disagreement (1/0)" = 1:2, "Prop. Disagreement" = 3:4),
       table = FALSE,
       include.proj.stats = FALSE,
       file = "tables/h2_langs_lpm_panel.tex")








